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ABSTRACT 

. Context. A hypothesis for sunspot formation is the buoyant emergence of magnetic flux tubes created by the strong radial shear at 

the tachocline. In this scenario, the magnetic field has to exceed a threshold value before it becomes buoyant and emerges through the 
whole convection zone. 

\ Aims. We follow the evolution of a random seed magnetic field with the aim of study under what conditions it is possible to excite 

I^T I . the dynamo instability and whether the dynamo generated magnetic field becomes buoyantly unstable and emerges to the surface as 

expected in the flux-tube context. 

Methods. We perform numerical simulations of compressible turbulent convection that include a vertical shear layer Like the solar 
tachocline, the shear is located at the interface between convective and stable layers. 

Results. We find that shear and convection are able to amplify the initial magnetic field and form large-scale elongated magnetic 
structures. The magnetic field strength depends on several parameters such as the shear amplitude, the thickness and location of the 
shear layer, and the magnetic Reynolds number (Rm). Models with deeper and thicker tachoclines allow longer storage and are more 
favorable for generating a mean magnetic field. Models with higher Rm grow faster but saturate at slightly lower levels. Whenever 
^ . the toroidal magnetic field reaches amplitudes greater a threshold value which is close to the equipartition value, it becomes buoyant 

■ and rises into the convection zone where it expands and forms mushroom shape structures. Some events of emergence, i.e. those with 

I ' the largest amplitudes of the initial field, are able to reach the very uppermost layers of the domain. These episodes are able to modify 

O \ the convective pattern forming either broader convection cells or convective eddies elongated in the direction of the field. However, 

^ . in none of these events the field preserves its initial structure. The back-reaction of the magnetic field on the fluid is also observed in 

(/^ ■ lower values of the turbulent velocity and in perturbations of approximately three per cent on the shear profile. 

' Conclusions. The results indicate that buoyancy is a common phenomena when the magnetic field is amplified through dynamo 

action in a narrow layer It is, however, very hard for the field to rise up to the surface without losing its initial coherence. 
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1. Introduction models, a set of complications have questioned the feasibility of 



CO 



this scenario. These may be summarized as follows: 



X 



Sunspots appear at the solar surface following a 11 -year cycle. 
CS| They reveal the presence of strong magnetic fields in the solar 1 . Magnetic flux tubes at the base of the convection zone should 
O interior, and suggest the existence of a dynamo process gov- have a strength of around 10^ G in order to become unsta- 
^ erning its evolution. However, the process by which sunspots ble a nd then cross the en tire convection zone up to the sur- 
^ are formed is still unknown. At the solar surface, sunspots are face (ICaUgari et al.l ll995) . This has been considered a prob- 
observed as bipolar patches of radial magnetic field. This intu- lem since the magnetic energy density corresponding to such 
itively suggest that they are formed by the emergence of horizon- field is one to two orders of magnitude larger than the kinetic 
tal concentrations of magnetic field lines, often called magnetic energy of the turbulent motions (the equipartition energy). 
^ flux tub es. Since P arker's original proposal of magnetic buoy- Note, however, that in the presence of shear this is not any 
ancy (Parked ll95"5l) . the model has evolved through the thin flux longer an upper limit for the amplitude of the magnetic field 
tube approximation (Spruit 1981), to numerical simulations of (see e.g. Kapyla & Brandenburg 2009). 
the emergence of 3D magnetic flux tubes (lFanll2008l) . The buoy- 2. Sunspots are coherent magnetic structures which means that 
ant rise of flux tubes from the tachocline up to the surface is the flux tubes should preserve their integrity while they rise 
currently the most widely accepted mechanism of sunspot for- through the entire convection zone. However, this region is 
mation. During the last three decades much work has been done highly turbulent and stratified, spanning more than 20 pres- 
in order to understand the buoyancy phenomena and to reconcile sure scale heights, so in addition to the 10^ G strength, the 
the results of flux tube models with phenomenological sunspots tubes must have a certain amou nt of twist in order to resist 
rules such as the Joy's law or the topological difference between the turbulent diffusion (see e.g. Emonet & Moreno-InsertisI 
the two spots in a pair. In spite of the fact that these basic ob- The current results do not conclu- 
servations have been qualitatively reproduced by the flux tube sively yield the amount of twist required for the tubes 
to rise coherently up to the surface. On one hand, 2D 
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vorticity generated by the buoya n t rise (ISchusslerl 



1979 



1996t 



Moreno-Insertis & Emonet 1996; Lon scope et alJ 
Emonet & Moreno-Insertis 1998; Fan et al. 1998^ On the 
other hand, rising 3D flux tubes require less twist thanks to 
the t ension forc es due to the longitudinal curvature of the 
tube (lFanll2QQlh . Three dimensional simulations in spheri- 
cal coordinates require fine tuning of the initial twist for the 
tubes to emerge to the surface with the observed tilt (iFanI 
2008). 

3. MHD simulations of convection in Cartesian coordinates 
have been a ble to produce larg e-scale magne tic fields 
through o? ( Kapvla et aP l2009bh and olVl ( Kapvl a et aP 
I2OO8) dynamo action. However, the magnetic fields gen- 
erated on these simulations are more homogeneously dis- 
tributed in space rather than in the form of isolated mag- 
netic structures. To the date, however, no direct simulations 
have been able to spontaneously form sunspot-like magnetic 
structures. 

4. Furthermore, instead of showing signs of buoyant rise and 
emergence of the magnetic field, the 3D dynamo simulations 
above have shown that the magnetic field is pumped down 
by convective downflows and tends to remain in the stable 
layer ( see also iTobias etaP 1 19981 12OOII: lOssendrijver et al.l 
l2002l) . lFanet al.1 (l2003l) studied the evolution of an isolated 
flux tube in a turbulent convection zone. They found that co- 
herent rise is possible as far as the magnetic buoyant force 
overcomes the hydrodynamic force from convection, i.e., 
obeying the condition > (i^p/a)^/^5eq, where Bq is 
the initial magnetic field strength. Hp is the pressure scale 
height, and a is the radius of the tube. In their case where 
(i^p/a)^/^ ^ 3, magnetic flux tubes with Bq > 35eq are 
able to rise coherently without being affected by convection. 
Similar results were found in a similar setup in spherical ge- 
ometry ( Jouve & Brun 2009). However, this does not address 
the issue with the pumping since, firstly, a strong magnetic 
flux tube is imposed on the convective layer, and secondly, 
magnetic pumping is an effect related with the g radient of 
turbulence intensity (Kich atinov & Riidigeiil 19920 . which is 
not present in these simulations. 



In addition to the thin flux tube approximation and simu- 
lations of rising flux tubes on stratified atmospheres, other re- 
cent attempts have been made in order to simulate the forma- 
tion of a magnetic layer thr ough the interaction of a n imposed 
she ar in convectiyely sta ble (IVasil & Brummelll2008b and unsta- 
ble (ISilvers et al.ll2009 a^ atmospheres, with an imposed vertical 
magnetic field. They have found that, unlike in the cases of im- 
posed toroidal magnetic layers, buoyancy instability is harder 
to excite whe n the magnet i c field is generated by the shear. 
More recently [Silvers et aP (l2009bl) . have found, with a sim- 
ilar setup, that the buoyancy may be favored by the presence 
of double-diffusive instabilities (these in turn depend on the ra- 
tio between the thermal and magnetic diffusivities, x/v^ often 
kn own as the inver s e Rob erts number). Later independent study 
of IChatterjee et al.l (l2010h has confirmed this result. However, 
in most of the current models, the presence of stratified turbu- 
lence, self-consistent generation of the magnetic field, or both, 
are omitted. 

In view of the above mentioned issues, other mechanisms 
have been proposed in order to explain sunspots. These are 
related to instabilities due to the presence of a diffuse large- 
scale magnetic field in a highly str atified turbulent medium 
( Kleeorin & Rosachevskii 1994; iRogachevskii & KleeorinI 
i2007.: .Brandenburg et al.. .2010a.b) . Mean-field models using 



this mechanism are able to produce strong flux concentra- 
tions but this has not yet been achieved in direct numerical 
simulations. 

Here we present numerical simulations of compressible tur- 
bulent convection with an imposed radial shear flow located 
in a sub-adiabatic layer beneath the convective region. For nu- 
merical reasons, explained below, we do not include rotation 
in our setup, i.e. the turbulence is not helical, and an aQ dy- 
namo is not expected. Nevertheless, recent studies have shown 
that mean-field dynamo action is possible due to non-helical tu r- 
bulence and shear (Brandenbu rg 20()5al lYousef et"aP 12008 Jbl : 
iBrandenburg et al.ll2008b . The nature of this dynamo is not yet 
entirely clear and may be attri buted to the so called shear- 
current effect (iRogachevskii & Kleeorin 2003, 2004) o r to th e 
incoherent, stochastic, a-effect ( Vishniac & Brande nburg 199"7l). 
According to ( Brandenburg et al..2008) the latter explanation is 
consistent with the turbulent transport coefficients. 

Based on these results, we expect the development of a mean 
field magnetic field, i.e. dynamo action, with a system that mim- 
ics, as far as possible, the conditions of the solar interior, specifi- 
cally in the lower part of the convection zone and the tachocline. 
As the shear is localized in a very narrow layer, we also expect 
the formation of a magnetic layer and the subsequent buoyancy 
of the magnetic fields. 

A similar setup was studied recently bylTobias et al.| (|2008h . 
They reported the appearance of elongated stripes of magnetic 
field in the direction of the shear. However, since they considered 
the Boussinesq approximation in their simulations, no buoyancy 
was observed. They also do not report the presence of a large 
scale dynamo. 

Two important features distinguish the simulations presented 
here from previous studies in the context of flux tube formation 
and emergence. Firstly, we consider a highly stratified domain 
with ^ 8 scale heights in pressure and ^ 6 scale heights in den- 
sity. Secondly, we do not impose a background radial magnetic 
field but allow the self-consistent development of the field from 
a initial random seed. 

Even though this is a complicated setup where it is difficult 
to analyze the different processes occurring independently, we 
believe that these simulations may give us some light on the 
current paradigm of sunspot formation. There are several im- 
portant issues that we want to address with the following sim- 
ulations. (1) What are requirements for dynamo action in the 
present setup? This includes the dependence of the dynamo ex- 
citation on several parameters such as the amplitude of the shear, 
thickness and location of the shear layer, and the aspect ratio of 
the box. (2) What is the resulting configuration of the magnetic 
field? In particular, whether the field is predominantly in small 
or large scales, and whether it is organized in the form of a mag- 
netic layer or isolated m agnetic flux tubes. (3) Is the buoyancy 
instability (iParkeJ 19551 ) operating on these magnetic structures? 
If yes, (4) how it depends on the parameters listed above? (5) Is 
it possible to have magnetic structures strong enough to emerge 
from the shear layer to the surface without being affected by the 
turbulent convective motions? (6) Finally, it is important to study 
how these strong structures back-react on the fluid motions, in- 
cluding the shear profile as well as the convective pattern. 

Another important issue that may be addressed in this con- 
text is the mechanism that triggers the dynamo instability. With 
the recent developments on the test-field method, it is possibly 
to compute the dynamo transport coefficients and have a better 
understanding on the underlying mechanism. 
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We have organized this paper as follow: in Sect.|2]we provide 
the details of the numerical model, in Sect. [3] we describe our 
results. We summarize and conclude in Sect. ID 

2. The model 



Our model setup is simila r to that used by e.g. 'Brandenburg et alj 
( 199 6) and Kapyla et al.l 42008). A rectangular portion of a star 
is modeled by a box whose dimensions are {L^^ Ly^ Lz) = 
(4, 4, 2)(i, where d is the depth of the convectively unstable layer, 
which is also used as the unit of length. The box is divided 
into three layers, an upper cooling layer, a convectively unsta- 
ble layer, and a stable overshoot layer (see below). The follow- 
ing set of equations for compressible magnetohydrodynamics is 
being solved: 

dA 

U X B — fJLQrjJ, 



dt 



DU 

~Dt 



--Vp- 
P 

Dt 



Dlnp 
Dt 

1 



1, 



JxB^-V -2iypS- 



P 



V ■ KVT ^2vS^ 



u-u 



(0) 



(1) 



(2) 



(3) 



(4) 



where D/Dt = d/dt -\- U - V is the total time derivative. A is 
the magnetic vector potential, S = V x A is the magnetic field, 
and J = V X S//io is the current density, /io is the magnetic 
permeability, r] and u are the magnetic diffusivity and kinematic 
viscosity, respectively, T is the temperature, s is the specific en- 
tropy, K is the heat conductivity, p is the density, U is the ve- 
locity, and g = —^z is the gravitational acceleration. The fluid 
obeys an ideal gas law p = pe(7 — 1), where p and e are the pres- 
sure and internal energy, respectively, and 7 = cp/cy = 5/3 is 
the ratio of specific heats at constant pressure and volume, re- 
spectively. The specific internal energy per unit mass is related 
to the temperature via e = cyT. The rate of strain tensor S is 
given by 



(5) 



where the commas denote derivatives. The last term of Eq. ^ 
relaxes the horizontally averaged mean velocity U towards a tar- 
get profile U^^^ where rf = 2y^dg) is a relaxation time scale. 
In a non-shearing case this value corresponds to ^ l/2rturn, 
where rturn = (^rms^f)"^ gives an estimate of the convective 
turnover time of the turnover time. The target profile is chosen 
so that it mimics the radial shear in the solar tachocline. More 
details of the shear profiles are given in Sect. (12.3) 1. 

The last term of Eq. describes cooling at the top of the 
domain according to 

T-T4 



■ cool 



(6) 



where f{z) is a profile function equal to unity in z > zs and 
smoothly connecting to zero below, and Fq is a cooling lumi- 
nosity chosen so that the sound speed in the uppermost layer 
relaxes toward T4 = T{z = Z4). 

The positions of the bottom of the box, bottom and top of 
the convectively unstable layer, and the top of the box, respec- 
tively, are given by (^1, ^2, ^3, ^4) = (—0.8, 0, 1, 1.2)d. Initially 
the stratification is piecewise polytropic with polytropic indices 
{mi^m2^ms) = (6, 1, 1), which leads to a convectively unsta- 
ble layer above a stable layer at the bottom of the domain and an 
isothermal cooling layer at the top. 
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Fig. 1. Vertical profiles of mean density, pressure, temperature, 
specific entropy and azimuthal velocity (Uy) in the initial state 
(solid lines) and the in the thermally saturated state (dotted). The 
dotted vertical lines at z = and z = d denote the bottom 
and top of the convectively unstable layer, respectively. The last 
two panels show the turbulent rms-velocity and the energy of the 
toroidal mean magnetic field in the saturated phase. 



2. 1. Nondimensional units and parameters 

Dimensionless quantities are obtained by setting 



(7) 
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where po is the initial density at ^2. The units of length, time, 
velocity, density, entropy, and magnetic field are 



[x]=d, [t]-- 
[s] = cp , [B] = \/~dgp^ . 



(8) 



We define the fluid and magnetic Prandtl numbers and the 
Rayleigh number as 



Pr=^ 

Xo 



Pm = - 



Ra = 



(9) 



where xo = K/ (pmCp) is the thermal diffusivity, and pm is the 
density in the middle of the unstable layer. For the magnetic dif- 
fusivity we consider a z-dependent profile which gives an or- 
der of magnitude smaller value in the radiative layer than in the 
convection zone through a smooth transition. Thus the magnetic 
Prandtl number in the stable layer is Pm =15 and in the con- 
vection zone Pm = 1.5. 

The entropy gradient, measured in the middle of the con- 
vectively unstable layer in the initial non-convecting hydrostatic 
state, is given by 



1 ds 
Cp dz 



V- V; 



ad 



(10) 



where V — Vad is the superadiabatic temperature gradient with 
Vad = 1-1/7, V = (SlnT/ainp),^,wherez^ = ^(^3+^2). 
The amount of stratification is determined by the parameter ^0 = 
(7 — l)cvT4/ {gd), which is the pressure scale height at the top 
of the domain normalized by the depth of the unstable layer. We 
use in all cases ^0 = 0.12, which results in a density contrast of 
about 120. We define the fluid and magnetic Reynolds numbers 
via 



Re 



Km : 



PmRe, 



(11) 



where u^ms is the rms value of the velocity fluctuations and /cf = 
27r/(i is assumed as a reasonable estimate for the wavenumber 
of the energy-carrying eddies. Our definitions of the Reynolds 
numbers are smaller than the usually adopted ones by a factor of 
2tt. The amount of shear is quantified by 



Sh: 



/7o/4 

^rms^f 



(12) 



where /7o is the amplitude and d^ the width of the imposed shear 
profile (see below). The equipartition magnetic field is defined 
by 

Beq = ( WC/'))^/^ (13) 

where the angular brackets denote horizontal average. For a bet- 
ter comparison between the magnetic and kinetic energies, we 
evaluate 5eq at the center of the shear layer, z = z^ef- 

The simulations were performed with the PENCIL COD^], 
which uses sixth-order explicit finite differences in space and 
third order accurate time stepping method. 



2.2. Boundary conditions 

In the horizontal x and y directions we use periodic boundary 
conditions and at the vertical (z) boundaries we use stress-free 
boundary conditions for the velocity. 



(14) 



For the magnetic field vertical field condition is used on the up- 
per boundary whereas perfect conductor conditions are used at 
the lower boundary, i.e. 



Bx = By 



B 



y,z 



B. 



0, {z = z^) 
0, {z = zi) 



(15) 
(16) 



respectively. The upper boundary thus allows magnetic helicity 
flux whereas at the lower boundary does not. This is likely to 
be representative of the situation in a real star where magnetic 
helicity can escape via the surface but does not penetrate into 
the core. 

2.3. Shear at the base of the convection zone 

In order to mimic the tachocline at the base of the solar convec- 
tion zone we introduce a shear profile 



^ref 



4 



(17) 
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where Zref is the reference position of the shear layer. Given the 
uncertainties of the radial position and width of the tachocline 
we perform parameter studies where Zref and d^ are varied. 
Furthermore, it is of general interest to study how dynamo ex- 
citation is depends on varying Uq and the ratio ofUo /d^. 

3. Results 

With the purpose of addressing the questions raised in the in- 
troduction we perform a series of simulations with the model 
described above where some properties of the shear layer are 
varied. We first study the hydrodynamic properties of the sys- 
tem and the conditions for dynamo excitation. The results of this 
parameter study are summarized in Table [T] Then we study the 
topological and buoyant properties of the magnetic fields gener- 
ated in some characteristic runs and compare them with models 
with different aspect ratio and higher resolution (see Table O. 
We finalize this section with the study of the magnetic feedback 
on the plasma motion and the computation of the turbulent coef- 
ficients that govern the evolution of large-scale magnetic fields 
in our simulations. 



3.1. Hydrodynamic instabilities 

In our simulation setup, two kinds of hydrodynamical instabili- 
ties may develop, namely the Kelvin-Helmholtz (KH) instability 
due to the imposed shear and stratification, and the convective 
instability due to the superadiabatic stratification in the middle 
layer. From hydrodynamical runs, we find that the convective 
instability develops early, at tUrmsh ~ 40, whereas the KH- 
instability develops at tUrmsh ^ 100 in runs with the strongest 
shear. After a few hundred time units the velocity reaches a sta- 
tistically steady state (constant rms- velocity), and the toroidal 
velocity achieves the desired shear profile. A thermally relaxed 
state (constant thermal energy), however, is reached only after a 
few thousand time units. This time depends on the radiative con- 
ductivity (K) but also on the imposed shear, which produces vis- 
cous heating that modifies the thermal stratification of the system 
as it may be seen in the top panels of Fig. [T] The final velocity 
profile which includes shear, convection, and the KH instability, 
does not allow the possibility of including rotation in the model. 
If it is done, the system develops mean field motions in the hor- 
izontal direction which are undesirable in the present study. 
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Table 1. Summary of the runs. 



Run 








Sh 


Re 


Rm 


Ma 


A[10-'] 


Beq 




By 


B/Bt 


SOO 


0.00 


0.05 


-0.05 


0.0 


8.9 


14.7 


0.028 


-2.16 


0.06 


— 




— 


SOI 


0.06 


0.05 


-0.05 


3.7 


8.2 


13.6 


0.026 


-1.29 


0.05 


— 


— 


— 


S02 


0.09 


0.05 


-0.05 


6.8 


13.4 


22.3 


0.042 


-0.34 


0.21 


— 


— 


— 


SOS 


0.12 


0.05 


-0.05 


7.3 


15.3 


25.5 


0.053 


0.03 


0.24 


— 


— 


— 


S04 


0.15 


0.05 


-0.05 


7.6 


19.9 


33.2 


0.062 


1.72 


0.40 


0.31 


5.32 


0.34 


DOl 


0.15 


0.05 


-0.05 


7.6 


19.9 


33.2 


0.062 


1.72 


0.40 


0.31 


5.32 


0.34 


D02 


0.15 


0.05 


-0.09 


7.7 


19.7 


32.8 


0.062 


1.72 


0.39 


0.40 


4.92 


0.36 


DOS 


0.15 


0.05 


-0.13 


7.6 


19.9 


33.1 


0.062 


1.36 


0.39 


0.42 


4.83 


0.36 


D04 


0.15 


0.05 


-0.17 


7.3 


20.1 


34.9 


0.066 


1.37 


0.39 


0.56 


6.15 


0.37 


TOl 


0.15 


0.05 


-0.05 


7.6 


19.9 


33.2 


0.062 


1.72 


0.40 


0.31 


5.32 


0.34 


T02 


0.15 


0.07 


-0.05 


6.0 


18.1 


30.1 


0.057 


1.00 


0.33 


0.34 


4.79 


0.36 


TOS 


0.15 


0.09 


-0.05 


4.9 


17.2 


28.7 


0.054 


0.91 


0.30 


0.31 


5.40 


0.37 


T04 


0.15 


0.11 


-0.05 


4.4 


15.8 


26.3 


0.049 


0.88 


0.26 


0.44 


6.15 


0.41 


T04o 


0.15 


0.11 


-0.05 


4.1 


16.6 


20.8 


0.052 


-4.3 


0.28 








T04t 


0.15 


0.11 


-0.05 


3.8 


18.1 


30.2 


0.056 


1.18 


0.36 


0.40 


3.75 


0.39 


TOS 


0.15 


0.15 


-0.05 


3.4 


15.0 


24.9 


0.045 


0.39 


0.24 








T06 


0.15 


0.20 


-0.05 


2.9 


13.1 


22.0 


0.041 


-0.04 


0.19 









Notes. Here Ma = Urms/igdy^^, A is the growth rate in units of rtum, ^rms = Brms/Beq, By = mdix{\By\) / Beq, where the equipartition 
magnetic field is evaluated at 2; = Zrei (see Eq.fTs]) and is given in units of y/dgpo/Io. Finally B/Bt = ^ {B)l/{B^) corresponds to fraction 
of mean magnetic field compared with the total field. The numbers are given for the saturated state of the dynamo. In all these simulations the 
resolution is 128^ grid points, Pr = 20 and Ra ^ 8.3 x 10^. 

Table 2. Summary of runs with different aspect ratio (upper two rows) and higher resolution (lower rows). 



Run 


{Lx, Ly, Lz) 


Uo/{dgy/' 


dz 




Pr 


Sh 


Re 


Rm 


Mak 


Mas 


A[10-'] 


Beq 


Brms 


By 


B/Bt 


AROl 


(8,4,2)d 


0.15 


0.11 


-0.05 


20 


4.4 


15.8 


26.3 


0.050 


0.050 


0.89 


0.25 


0.50 


6.08 


0.40 


AR02 


(4,8,2)(i 


0.15 


0.11 


-0.05 


20 


3.8 


18.4 


30.6 


0.058 


0.058 


0.30 


0.32 


0.22 


4.73 


0.29 


DOSb 


(4,4,2)d 


0.15 


0.05 


-0.10 


12 


7.9 


31.8 


47.7 


0.062 


0.060 


2.57 


0.43 


0.31 


4.16 


0.31 


T04a 


(4,4,2)0? 


0.15 


0.1 


-0.05 


20 


5.1 


25.0 


37.6 


0.047 


0.047 


2.23 


0.26 


0.41 


5.06 


0.36 


T04b 


(4,4,2)(i 


0.15 


0.1 


-0.05 


12 


4.9 


25.6 


38.5 


0.052 


0.048 


2.19 


0.27 


0.43 


4.80 


0.37 


T04c 


(4,4,2)(i 


0.15 


0.1 


-0.05 


3 


4.4 


28.6 


43.0 


0.056 


0.054 


1.21 


0.30 


0.35 


4.13 


0.36 



Notes. Most of the quantities here are defined in Table [T] We have added here the scales of the domain, the Prandtl number and the Mach number 
at the kinematic (Mak) and saturated (Mas) stages. Models with different aspect ratio have the same spatial resolution than models in Table [T] In 
Runs DOSb and T04a-c the resolution is 256^ grid points. 



3.2. Dynamo excitation 

In order to study the dynamo excitation we follow the evolution 
of an initial random seed magnetic field of the order of 10~^5eq- 
Since there is no rotation _and becausethe vertical component of 
the large-scale vorticity Wi = CijkUj^k ^ 0, the average ki- 
netic helicity and mean-field a-effect are expected to vanish (e.g. 
iKrause & Radleilll98 0). However, it is still possible to excite 
a large-scale dynamo due to non-helical tu rbulence and shear 
(lBrandenburgll2QQ5al; lYousef et al.ll2QQ8allbl: iBrandenburg et aP 
'2008). In this part of this study we vary three parameters: the dif- 
ferential rotation amplitude, Uq (Runs S00-S04), the location of 
the shear layer in the convectively stable region, z^ef (Runs D02- 
D04), and the thickness of the shear region, dz (Runs T02-T04). 

We note that no small-scale dynamo is excited in the non- 
shearing case. With the runs in Set S we find that the critical 
value of shear above which a dynamo is excited is Sh 7 (see 
the top panel of Fig. O. The value of Sh is roughly three times 
larger than that computed with typical values of the rms-velocity 
and length scales in the solar tachocline0. 

^ We have considered AQ — 33nHz, rtac = 0.7/ ^0, dtac = 
O.O5i^0 (e.g. IChristensen-Dalsgaard & ThompsonI I2007L and refer- 
ences therein), i^rms = 10^ cm s~^, and = 27r//, with / = 10^° 



In Set D, we move the shear layer deeper down into the sta- 
ble layer, where the turbulent diffusivity is expected to be less, 
thus allowing a longer storage of the magnetic fielcfl. From the 
hydrodynamical point of view, these simulations have a convec- 
tion zone that, due to viscous heating, extends into the initially 
stable layer in the saturated state as it can be seen in the ver- 
tical profile of entropy corresponding to Run D04 (red line) in 
Fig. [T]). The three bottom panels of the same figure (see also 
Table [B show (from topjo bottom) the vertical profiles of the 
mean toroidal velocity (Uy), the turbulent rms-velocity and By 
in the relaxed state. Through this paper we consider averages 
taken first over the y direction and then the averaging over the 
other directions is performed. Note that the iirms -profiles have 
been computed neglecting the mean flow. The curves of the tur- 
bulent velocity in the models that include the shear layer exhibit 
a bump where the velocity gradient is located. This feature does 
not appear in non-shearing models (see dotted line), suggesting 
that turbulent motions are developing in this regions. The bump 
is more pronunced in models with larger Sh which indicates that 



cm (see Table 1 of IBrandenburg & Subramanianll2005h . which result in 

Shtac ~ 2. 

^ Note that Run S04 is the same as Runs DOl and TOl. 
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Fig. 2. Time evolution of the rms value of the total magnetic field 
normalized with the equipartition value at the center of the shear 
layer. Different lines/colors correspond to different runs as indi- 
cated in the legends in each panel. 

this turbulent motions are probably due to the KH instability. 
This difference in iirms reflects in the value of the equipartition 
magnetic field (Eq.[T3l) which differs from model to model. 

We find that the amplitude of the magnetic energy increases 
for deeper tachoclines. This is an expected result since below 
the convection zone the magnetic diffusivity has smaller values 
which allows a larger storage. The vertical distribution of the 

2 

toroidal component of the magnetic energy (5 ) peaks roughly 



at the center of the shear layer. The tail of this curve towards 
the convection zone may hint how buoyant the magnetic field 
is on each model. We expect that the larger the magnetic field 
(Run D04) the more magnetic flux becomes buoyant and rises 
into the upper layers. 

The time evolution of the magnetic field of the runs in Set D 
are shown in the middle panel of Fig.O We find that the growth 
rate of the magnetic field, A = dln^rms/dt, decreases slightly 
as the shear layer is moved deeper. This is not very clear in 
Table [T] since the error of this quantity is of the same order as 
the difference between different runs. However, it is clear from 
both, the figure and the table, that the volume averaged magnetic 
field increases with the tachocline depth. 

In the third group of simulations (Runs T01-T04), we in- 
crease the width of the shear layer gradually from = 0.05(i in 
Run TOl to (i^ = 0.11(i in Run T04. This change implies lower 
values of the shear parameter Sh (Eq.[T2]), but also larger fraction 
of the tachocline in the turbulent and stable layers. The results, 
presented in Table[T]and depicted in Figs. [T] and [2] (see legends), 
indicate that smaller values of the shear, Sh ^ 3, are still able to 
excite the dynamo at the price of a lower growth rate. We notice 
that the growth rate depends on the magnetic Reynolds number 
(compare Runs T04 and T04b). The fact that a smaller shear 
generates a larger magnetic field is a counter-intuitive result and 
does not agree with pre vious mean-field studies on the thicknes s 
of the solar tachocline dGuerrero & de Gouveia Dal Pinoll2QQ7l) . 
However, as it will be explained below, the important fact is that 
this configuration seems to be favorable to a longer storage of 
magnetic field in the stably stratified layer, with the advantage 
that here the effects of the shear on the thermal properties of 
the fluid are less important than in the previous sets of simula- 
tions. With the settings of Run T04 we find that the critical mag- 
netic Reynolds number is between 20.8 (Run T04o) and 26.3 
(Run T04o). 

Based on the results above, we perform another series of 
simulations whose parameters and results are summarized in 
Table O Runs AROl and AR02 correspond to Run T04 but with 
aspect ratios [L^^Ly^Lz) = (8,4,2)(i and (4,8,2)(i, respec- 
tively. Runs D03b, and T04a to T04c have essentially the same 
configuration than Runs DOS and T04 with 256^ grid points res- 
olution. The values of u and K have been modified in order to 
obtain different Reynolds and Prandtl numbers. 

The model with larger extent perpendicular to the direction 
of the shear does not show differences with respect to the ref- 
erence case. On the other hand, the model with larger extent 
in the direction of the shear results in a reduced growth rate 
and in smaller values of the volume averaged magnetic field 
(^rms), the maximum amplitude of the toroidal magnetic field 
and of the large-scale field. The reason for these changes is the 
increa se of the Urmp veloc ity, which diminishes the effective 
shear. lYousef et"aD (l2008ah have obtained that both the growth 
rate and the scale where the maximum magnetic energy is con- 
centrated, converge to a value which is independent of the rele- 
vant length scale (Lz in their case). Such convergence analysis is 
computationally very expensive to be performed for our system. 
The models with higher resolution (see Runs D03b and T04b in 
Table [2] and bottom panel of Fig. [2]) show a larger growth rate 
when compared with their corresponding low resolution cases. 
These runs, however, saturate at sligthly smaller amplitudes of 
^rms- Similar weak dependence o n Rm has be en reported from 
simulations with horizontal shear (iKapyla et a l. 2010a). 

We have also changed the Rm by considering a different in- 
put heat flux (i.e. different Pr). The results of Runs T04a-T04c 
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Fig. 3. Snapshots of the vertical velocity (top row) and the 
toroidal magnetic field {By) in different planes of the domain, 
for Runs T04c (from left to right) in the kinematic phase. The 
vertical slices correspond to yz (left side of the box) and xz 
(front of the box) planes. The top plane corresponds to the hor- 
izontal boundary between the convective and the cooling layers 
at z = d and the plane shown below the box corresponds to the 
base of the unstable layer (z = 0). 
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show that A depends on Pr. There is no considerable difference 
between the rms magnetic field of these runs. 

3.3. Morphology of the dynamo generated magnetic field 

In the initial stages of evolution, the magnetic field is dominated 
by small scales and although it is possible to distinguish the 
stretching effects of the shear, the structures formed are small 
compared with the size of the box. In Fig. [3] we present snap- 
shots of the Run T04c in the kinematic phase. 

Mean magnetic fields take longer to develop in all simula- 
tions but in the saturated state they correspond to a considerable 
fraction of the total field. This_can be seen in the upper panel 
of Fig. m where the ratio 6rms/ ^rms for Run S04 is shown as a 
function of time. It is important to notice that mean values refer 

here to averages over the y direction, i.e. 5rms = ^(^)^, and 

6rms = ^ {B'^)y ~ (^)r perform average over the 

X direction because the toroidal magnetic field varies also in x 
(see below). The average over x and time is performed after the 
mean and fluctuation values are computed. 

In this figure, the dashed line corresponds to the average 
value of this ratio (1.9) in the saturated phase of the dynamo, 
indicating that the mean magnetic field is roughly a third of the 
total magnetic field. For the sake of clarity only a single run is 
shown in this panel, but similar results are obtained for all simu- 
lations in Sets D and T. 

In the bottom panel of Fig.|4l the vertical profiles of the mean 
(solid lines) and fluctuating (dashed lines) magnetic fields for 
four representative runs are shown. We find that the mean mag- 



Fig. 4. Upper panel: ratio between small scale, 6rms, and mean, 
^rms, magnetic fields for Run S04. The dashed line corresponds 
to the mean value of this ratio in the saturated phase. Lower 
panel: radial distribution of the mean (thick lines) and fluctu- 
ating (thin lines) magnetic field for the Runs S04, D04 and T04 
in the saturated phase. 



netic field is mainly located in the shear region, whereas the tur- 
bulent field, on the other hand, is more spread out inside the 
convection zone. It is noteworthy that for a thick tachocline the 
mean field is almost comparable with the fluctuating one, cover- 
ing a fraction of the stable layer where the fluctuations are weak 
(see continuous blue line). Simulations with higher resolution 
(see blue lines with diamond symbols in Fig. [?]) exhibit fluc- 
tuating magnetic field with vertical distribution and amplitude 
almost identical to the lower resolution case. The vertical profile 
of the mean field is roughly the same as in the lower resolution 
case but of smaller amplitude. 

The structure of the magnetic fields depends strongly on the 
structure of convection. In Fig. [5] we show snapshots of verti- 
cal velocity (JJz) and azimuthal magnetic field {By) for arbitrary 
times in the saturated state for Runs S04, D04, and T04 (from 
left to right, respectively). From the bottom panels of this fig- 
ure, it is possible to distinguish that at the base of the convec- 
tion zone, the toroidal field organizes in elongated structures of 
both polarities which span all across the azimuthal direction if 
the penetrative downflows are less intense. These stripes coex- 
ist with more random magnetic fields in regions located where 
jet-like overshooting is able to reach the stable layer. The size of 
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the magnetic structures is at least equal or larger than the scale 
of the convective eddies. This is more evident for the stream wise 
direction where the field occupies the full extent of the domain. 
The strong toroidal fields are mainly confined in the shear re- 
gion. This can be seen in the bottom panels of Fig. [5] and also in 
the bottom panel of Fig. [H where it is clear that the curve cor- 
responding to Run T04 (blue line) has a broader profile than the 
one corresponding to Run S04 (black line). The presence of the 
KH instability is evident in the yz plane where wave like struc- 
tures are observed. If the tachocline is thicker, the effects of the 
KH-instability are weaker. Another advantage of thicker shear 
layers is that they do not affect the thermodynamical properties 
of the fluid as strongly as in Sets S and T. The disadvantage is 
that they produce broad and more diffuse magnetic fields which 
do not totally agree with the picture of a flux tube. This is dis- 
cussed in more detail in the next section. 

Th ese results are in contradiction with those of iTobias et al.l 
(l2QQ8h who used a similar setup in the Boussinesq approxima- 
tion. We believe that the lack of a large-scale dynamo in their 
simulations may be due either to the wider shear profile that 
they used or to the smaller amplitude of the shear parameter. 
However, a visual inspection of their figure 1 1 suggests that a 
large-scale dynamo could indeed exist in the horizontal plane. 

3.4. Buoyancy 

According to linear theory (ISpruit & van BallegooijenI [l982l) , 
thin magnetic flux tubes become buoyantly unstable when the 
condition {35 > — 1/7 is fulfilled. Here /3 = is the 

plasma beta, S = \/ — Vad is the superadiabaticity (see Eq.fTOl), 
and 7 is the ratio of specific heats. This condition is satisfied in 
the convection zone where S > 0. From this relation it follows 
that below the convection zone, where ^ < 0, thin magnetic flux 
tubes require much higher field strength to become unstable. For 
a magnetic layer, the necessary and sufficient condition for the 
development of two-dimensional i nterchange mod es, in which 
we are interested here, is given by (lNewcomblll96lh : 



We find that for the cases with thinnest shear layer, the instability 
region is above the center of the shear layer where the toroidal 
magnetic field reaches its maximum value (see Fig. [6]). This im- 
plies that a large fraction of the magnetic field there becomes 
buoyantly unstable very quickly. For the cases with a thick shear 
layer (Set T), the magnetic field is unstable only within the con- 
vection zone. The magnetic field distribution of the Run T04 in 
Fig. [4] (blue line), indicates that a large fraction of the magnetic 
field is buoyantly stable. This configuration may explain why 
in these cases larger mean magnetic fields develop with smaller 
shear parameters than in sets S and D are used. 

Note, however, that this conclusion is based on averaged 
quantities. The evolution of local magnetic structures is rather 
complex and possibly also affected by the KH instability. The 
radial velocity gradient in Run T04 is smaller, so it is expected 
to have a less efficient KH instability and this should allow a 
longer stay of the magnetic field in the stable layer. This is hard 
to demonstrate since it is difficult to disentangle these effects in 
the simulations. Nevertheless, it seems that the longer the mag- 
netic field stays at the stable region the great er the final mean 
magnetic field strength. Hence, according to "FanI (I2QQ1 ), the 
number of scale heights the magnetic flux concentrations may 
rise, and its final structure depends only on how strong the field 
is in comparison to the turbulent convective motions. 



In the simulations it is observed that when the mag- 
netic field rises it expands as a consequence of decreasing 
density. The morphology of the magnetic fields in the xz 
plane (see bottom panels of Fig. O is reminiscent of the 
mushroom- shape that has been obtained in several 2D (.Schussled 
1979; Moreno-Insertis & Emonet 1996; Longcop eet al.lll996t 
Fan et al. 1998; Emonet & Moreno-Insertis 1998) and 3D (iFanI 
2001; Fanetal. 2003; Fan 2008; Jouve & Brun 2009) simula- 
tions of flux tube emergence. Such expansion may result in the 
splitting and braking of the tube. We do not impose any twist 
on the tubes since the dynamo generates the magnetic field self- 
consistently and the possibilities span from events that do not 
rise at all to events where buoyant magnetic fields rise up to the 
very surface. In the horizontal yz plane we observe that the mag- 
netic field lines remain horizontal in the stable layer where they 
form. When the magnetic field rises, the field lines bend in the 
convection zone, with their rising part in the middle of upward 
flows and other parts withheld to the downward flows. 

In Tables [T] and [2l we list the maximum values of By in the 
simulations. In most of the cases we find m8ix{By) ^ 55eq, 
and only in the models with a thicker (Runs T04 and AROl) 
or a deeper shear layer (Run D04), ma x(^^) can be some- 
what greater than 65eq. iFan et al.l (12003 ) argue that if Bq > 
{Hp I (i)^I^Bq^, Bq being the field strength in the shear layer and 
Hp the local pressure scale height, a buoyant flux-tube of radius 
a may rise without experiencing the convective drag force. Here, 
although the shape of the magnetic field is not a tube, we com- 
pute the same condition using a = dz. We obtain values of Bq 
going from 3.45eq in Run D04 to 2.25eq in Run T04 (2.35eq in 
higher resolution runs). This indicates that the dynamo generated 
magnetic field is insufficient to reach the surface without being 
modified by the convection. There are, however, a few cases in 
Runs D04, T04 and AROl, where strongest magnetic fields are 
able to reach the surface. When this occurs the magnetic field re- 
acts back on the flow and modifies the convective pattern. Broad 
convection cells elongated in the y direction are formed. In Fig.|7] 
we present one of these events for Run AROl. In the upper lay- 
ers, the field lines show a turbulent pattern except in the center 
and the right edge of the box where the field lines at the surface 
cross the entire domain in the y direction. In the higher resolu- 
tion cases, we do not observe events where the field lines, in the 
uppermost layers, remain horizontal accross the whole toroidal 
direction. Nevertheless, the expanding magnetic field may reach 
the surface forming broad convection cells (at least two times 
larger than the regular cells). One of these cases, correspond- 
ing to Run T04c is shown in the bottom panels of Fig. [71 In 
this event two stripes of magnetic field of opposite polarity are 
reaching the surface simultaneously, the magnetic field lines of 
these ropes reconnect forming a large scale loop oriented in the 
x-direction. 

In order to compute the rise speed of the magnetic field in 
these cases we have used the horizontal averages of By and con- 
structed the time-depth "butterfly" diagram shown in Fig. [8] The 
toroidal field is amplified in the shear layer, below z = (see 
dotted line). When it becomes buoyant it travels through the con- 
vection zone. The tilt observed in the contours of By may give 
a rough estimate of the vertical velocity. In the same figure we 
have drawn white dashed lines to guide the eye. For the emer- 
gence events in Run D04 (AROl), the estimated rise speed is 
Ub ^ 0.034((i^)i/2 0.030((i^)^/2). These values are smaU in 
comparison to the Urms of the models (see Tables [T] and O which 
does not agree with the rise of a flux tube in a stratified atmo- 
sphere. For instance, iFan et al ] (l2QQ3») obtain a final rise veloc- 
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Fig. 5. Same as Fig. [3] but for Runs S04, D04 and T04 (from left to right) in the saturated phase. Movies corresponding to these 
runs may be found at: www.nordita.org/~guerrero/movies. 
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Fig. 6. Buoyancy instability condition, Eq. ([TSb evaluated for Runs S04, D04 and T04. Continuous (dashed) lines correspond to 
the left (right) hand side of the equation, respectively. Dotted lines correspond to the interfaces between convectively stable and 
unstable layers, and to the center of the tachocline in each simulation. 



ity ^ 5 times larger than the rms-velocity. We should point out 
that in the simulations here, the strong shear increases i^rms by 
a factor of two or even more. If we compare the rise velocity of 
the magnetic field with the rms-velocity of a non- shearing case 
(Run SOO), we find that 2x5 is slightly larger. Magnetic pump- 
ing effects, which are evaluated below, may also play a role in 
braking the magnetic flux and in defining its final velocity. 

The r esults mentioned above d o not fully agree with those 
found by IVasil & Brummelll (l2QQ8h . where a toroidal magnetic 
layer is generated from a purely vertical field through an im- 
posed radial shear layer. Although their stratification is smaller 
than in our case, the buoyancy instability is hard to excite in their 
simulations and the buoyancy events are slowed down rather 
quickly. The lack of a dynamo mechanism, able to sustain the 
magnetic field, is maybe the reason of its rapid diffusion. 

3.5. Back reaction 

Temporal analysis of helioseismic data has revealed that the so- 
lar angular velocity varies by around 5 per cent with respect to 
its mean value. This fluctuation pattern, like the sunspots, fol- 
lows an 1 1 -year cycle, which suggests that it corresponds to the 
back reaction of the magnetic field on the plasma motion (e.g. 



iBasu & Antiall2QQ3l:lHowe et al.ll2009 l). It is also known that the 
amplitude of the meridional circulation varies with the amplitude 
of the magnetic field, being smaller when the cycle reaches its 
maximum (see e.g Basu & Antia 2003). Apart from large scale 
effects, the magnetic field may also affect the local properties of 
the plasma. Th i s is po ssible the case in the downflows found by 
Hind man et al.l (l2009l) nearby the bipolar active regions. 

In the simulations presented here, due to the different nature 
of the dynamo and to the numerical setup, there are neither peri- 
odic oscillations of the magnetic field nor large-scale circulation. 
We are able, however, to distinguish the differences between the 
plasma properties in the kinematic and the saturated stages. We 
also follow the deviations of the averaged flow occurring during 
peaks of the magnetic field amplitude. 

We find that also the turbulent rms-velocity, as well as the 
imposed large-scale shear flow, are affected by the magnetic 
field. In the top panel of Fig. [9] we show the vertical profile of 
the horizontal and time averaged Urms at the kinematic (dotted 
line) and saturated (continuous line) regimes of the dynamo. A 
decrease in the amplitude is observed, especially in the region 
where the magnetic field is more concentrated. The volume av- 
eraged Urms also changcs accordingly by ~ 5 per cent, as indi- 
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Fig. 7. Vertical velocity (left panel) and magnetic field lines (right panel) for one of the emergence events in Runs AROl and T04c. 
The colors of the magnetic field correspond to the values of By/B^ci- 



cated in Tabled The dashed line in the top panel of Fig. [9] shows 
that it further decreases in the shear region. 

In the middle and bottom panels of Fig.[9]butterfly_diagrams 
of the toroidal velocity fluctuations, SUy = Uy — Uy, in the 
z-time plane are shown for the Runs AROl and T04b, respec- 
tively. The dashed white line indicates the center of the shear 
layer. Positive values of this quantity below (above) the dashed 
line indicate faster (slower) toroidal velocities with respect to 
its mean value. Accordingly, negative values below (above) the 
dashed line indicate deficiency (excess) of the local y flow. The 
continuous white line corresponds to the time evolution of Brms 
(for clarity normalized to its maximum value). We show this in 
order to demonstrate that the changes in Uy depend on the am- 
plitude of the magnetic field. The overall result of the contour 
plots is that during the maximum of the magnetic field the shear 
is slightly stronger. In these diagrams, as well as in all of the 
simulations listed in Tables [T] and [2] (except Run T04r), the de- 
viation around the averaged value of Uy is found to be ~ 3 per 
cent. 

The deviation of the mean shear profile depends on the 
forcing time scale Tf in Eq. [3] In Run T04r, we consider the 
same settings that in Run T04 but varying Tf from 2{dg)^^'^ to 
10((i^)^/^ (^ 1 to ^ 5 turnover times, respectively). The results 
show fluctuations of ^ 10 per cent of the shear velocity, and 
spread out now to a wider fraction of the domain (see bottom 
panel of Fig.O. Since magnetic buoyancy is present also at the 
time of the maximum, the toroidal velocity in the middle of the 



convection zone is also affected. This change is not so strong as 
the one in the shear layer but is also a general characteristic of 
the simulations. 

The perturbations shown in Fig. [9] are reminiscent of the 
equatorial branches of the solar torsional oscillations which 
have a leading excess in the zonal flow, followed by a slightly 
less pronounced deficit of it. Similar results have also been 
found in global num erical simulations of oscillating dynamos by 
iBrown et aP (l201ll ). They have found that the branches of zonal 
flows follow the evolution of a poleward dynamo wave. 



3.6. Turbulent transport coefficients 

In order to study the mechanism that generates a mean-field dy- 
namo in our simulations, we compute the turbulent transport 
coefficients from a representative simulation (Run T04). In or- 
der to do this, we use the test-field me thod^ ntrodu ced in the 
geodynamo context by ISchrinner etaD ^2005", "200?) and cur- 
rently available in the PE NCIL CODE (e.g. Brandenburgll20Q5bt 
iBrandenburg et al.ll2008b . 

In mean- field theory, the electromotive force S = 
u X B governs the evolution of the large-scale magnetic field 
(Krause & Radler 1980). Under the assumption that the mean 
field varies smoothly in space and time, and that there is no small 
scale dynamo action, the electromotive force may be written in 
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Fig. 8. Butterfly diagram in the z-time plane for the toroidal 
magnetic field, By. The dotted lines corresponds to the top 
(z = d) and base {z = 0) of the convection zone. The tilted 
dashed lines show approximate trajectories of buoyant magnetic 
fields. A rough estimate of rise speed is computed with these 
lines as indicated. The event at t = 9020iirms^/ in the upper 
panel of Fig. [8] is also shown in Fig. [71 



terms of the large-scale magnetic field: 



aijBj 



Vijk J k ' 



(19) 



Considering large-scale fields that depend only on z we need 
four independent test fields in order to compute the 4+4 compo- 
nents of an and m^k isee sl detailed description of the method 
in Brandenb urg et al II2OO8I) . The novel feature of the test field 
method is that the test fields do not act back on the flow 
and that the turbulent diffusivity can also be computed, thus 
avoiding many of th e problems that plaque other methods (cf . 
iKapvla et al.l l2QlQbh . 

We discuss our results in terms of the following quantities: 



Vt = ^(^11 - 

7 = ^(<^21 



ayy = 0^22, 



•^22), Cry = liVll -^22), 

- ai2), = 4(^21 + 0^12), 



(20) 
(21) 
(22) 



which represent the inductive (a), diffusive (r]t), and pump- 
ing (7) effects of turbulence, respectively. In Fig. [TO] we show 
the vertical profiles of the kinetic helicity (u • cj), where = 
V X n is the vorticity, and the turbulent transport coefficients 
of Eqs. ([2Qb-([22b. normalized with the first order smoothing ap- 
proximation (FOSA) quantities: 



ao 



Vto 



3^f • 



(23) 



For computing the coefficients in Fig.[TOlwe use Run T04 in the 
purely hydrodynamic state. The magnetic diffusivity considered 
for the test fields is one order of magnitude larger that that used 
for the magnetic field in the original run. However, the magnetic 
Reynolds number (Rm ^ 4) is still sufficiently large to yield 
reasonable results (iKapyla et al.ll2Q09aD . 

We find that certain amount of helicity is generated in the 
system, but that it is still statistically consistent with zero. The 
coefficient axx^ despite large fluctuations, has a negative sign. 
This component of a contributes to amplifying the y compo- 
nent of the magnetic field, but its contribution is likely negligible 
when compared with that of the shear in the present models. 




Run AEOl, 




Fig. 9. Top: time averaged vertical profile of u^ms in the kine- 
matic (dotted line) and dynamical (continuous lined) stages for 
the Runs T04b. The dashed lines show the deviation from the av- 
erage at a peak in the magnetic field amplitude. The three lower 
panels are butterfly diagrams of the angular velocity deviation, 
8Uy = Uy - Uy, foY thc Runs HR01,T04b and T04r. The con- 
tinuous white line in these panels indicate the time evolution of 
Brms normalized to its maximum value in the time interval. 



The component ayy shows a zero mean value, but its vari- 
ance has a large amplitude, especially in the middle of the c on- 
vection zone. According to lVishniac & Brandenburgl(ll997l) . an 
incoherent a effect, zero in average but with finite variance, may 
generate enough inductive effects to sustain the dynamo. This 
suggest that the incoherent a effect may be the mechanism sus- 
taining the dynamo. 

The shear-current effect, arising from the inhomogeneity 
of the turbulence and the mean shear flow, may also resu lt in 
the generation of a mean field (Rogachevskii & Kleeorin 2 0031 
2004). The existence of this effect has been studied with models 
of forced turbulence and latitudinal shear by (Brandenburg et al.l 
I2OO8I: iMitra et al.ll2009l) . but no conclusive evidence has been 
found. Similar results have also been found from convection 
with horizontal shear ( Kapyla et al. 2009a). In the case presented 
here, with vertical shear, the relevant coefficient for this effect 
would be 7^13, which is not captured by the current version of 
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Fig. 10. Normalized profiles of the kinetic helicity, and the tur- 
bulent transport coefficients axx^ (^yy, 7 and r]t (from top to bot- 
tom) computed with the test-field method. 



the test-field method with only z-dependent test fields. Thus, we 
leave the study of the shear-current effect in turbulent convection 
with vertical shear for a forthcoming work. 

The turbulent diffusivity, r]t, is roughly four times larger than 
the FOSA estimate. When compared with the molecular dif- 
fusivity, it is also four times larger (i.e. of the same order as 
Rm ^ r]t/r]). This result agrees with previo usly computed tur- 
bulent diffusivity for convection (iKapyla et al. 2009a). The ver- 
tical profile of rjt differs somewhat from previous results, indi- 
cating the action of the radial shear on the turbulent diffusion. 
Another inter esting p oint is that is not close to zero, as in 
iKapvla et al.1 (l2009ah . This is likely explained by the stronger 
shear used here. This is also an indication of the tensorial char- 
acter of r]t and suggests that the distinct components of the mag- 
netic field may diffuse differently. 

Finally, the vertical profile of the turbulent pumping (7) in 
Fig. [TOl shows that there is downwards transport of the mean 
magnetic field at the bottom of the convection zone and upward 
transport close to the top. The pumping velocity is comparable 
with Urms but it is not enough to retain the magnetic field within 



the stable layer, which indicates that buoyancy is a more effi- 
cient mechanism than magnetic pumping in transporting mag- 
netic fields. Note, however, that the amplitude of this effect de- 
pends on the gradie nt of the level of turbulence, 7 = — Vt^t, 
where tjt = r] -\- r]t (* Kichatinov & Rtidigedl 19921) , which is not 
very large in the current simulations. 



4. Conclusions 

We perform numerical simulations of turbulent convection with 
a thin vertical shear layer with the aim of mimicing the condi- 
tions at the solar convection zone and the tachocline. This layer 
is located below the interface between a convectively stable and 
unstable layers. We consider the velocity gradient, the depth and 
the thickness of this layer as free parameters of the model and 
concentrated on two phenomena relevant for the solar cycle: dy- 
namo process and magnetic buoyancy. 

We find that it is possible to excite a large scale dynamo by 
the combined action of turbulent convection and a localized ra- 
dial shear. The efficiency of the dynamo (i.e., the growth rate), 
and the maximum amplitude of the magnetic field depend on the 
shear amplitude, the thickness of the shear layer, the aspect ratio 
and the magnetic Reynolds number of the system (see Fig.[2]and 
TablesIIlandO. 

The magnetic field is organized in elongated structures in 
the direction of the shear in the shear layer below the base of the 
convection zone. The large-scale magnetic field can have either 
polarity and may be up to the 40 per cent of the total field (solid 
lines in Fig. (4]). It coexists with turbulent magnetic field which 
is distributed all across the convection zone. In the stable layer 
the magnetic field exists in the regions of strong downflows. For 
fixed depth and thickness (Set S) a smaller growth rate and larger 
B^ms is obtained for a larger shear parameter. A critical value of 
Sh ^ 7 is required to excite the dynamo. 

If the tachocline is located deeper (simulations in set D), the 
magnetic field develops mainly in the stable layer which allows 
a longer field storage. In these cases the dynamo grows faster 
and with a larger fraction of mean field magnetic field. This 
configuration, however, is not optimal since with the numeri- 
cal resolution used here the thermodynamical properties of the 
fluid are affected by viscous heating in the shear region. If the 
depth of the shear layer remains constant but the thickness in- 
creases (Set T), the magnetic field grows slowly but it contains 
a larger fraction of mean magnetic field. This may be a conse- 
quence of the smoother shear profile which renders the magnetic 
field buoyantly unstable only in the superadiabatic part of the do- 
main (see right panel of Fig. O. A lower critical shear number is 
found to be able of excite the dynamo instability in these cases 
(Sh ^ 3.5). 

Larger magnetic Reynolds numbers are achieved in two 
ways, first by doubling the initial resolution and second by 
changing the input flux (i.e., different Pr number). On one side, 
for a fixed Pr and larger Rm, the dynamo exhibits a much faster 
growth. On the other hand, for a fixed resolution and different 
input flux, the growth rate changes proportionally to Pr. In these 
higher resolution Runs the final value of Brms is slightly be- 
low that in the corresponding cases with lower resolution. This 
may correspond to the dependence of the saturation value of the 
magnetic field with the magnetic Reynolds number or may also 
be the result of insufficient statistics. 

Since the system is not rotating, there is no kinetic helicity 
and hence no a-effect. The test-field analysis suggests that the 
probable mechanism triggering the amplification of the magnetic 
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field could be incoherent a-shear dynamo. However, other pos- 
sibilities like the shear-current dynamo can not be ruled out for 
the time being. 

Magnetic buoyancy is observed in all of the simulations. 
Based on horizontal averages we analyze the instability condi- 
tion for 2D interchange modes (Eq.[T8]and Fig. O. We find that 
in models with a deeper tachocline the buoyancy instability de- 
velops even in the stable layer whereas in models with thicker 
shear layers the magnetic field is unstable in the convection zone 
only. This, however, does not mean that there are not emergence 
events. As far as the toroidal magnetic field is strong enough it 
rises throw the convection zone forming mushrooms shape struc- 
tures. 

When the buoyant magnetic field is weak it is strongly mod- 
ified by local convective flows (see Fig. O. On the other hand, 
magnetic fields of strong amplitude formed at the shear layer are 
able to rise up to the surface and modify the convective pattern 
(Fig. [7]). They may form either large convection cells or convec- 
tive rolls may occupy all of the ^-extent. These clearest events 
are observed when the magnetic field in the shear layer exceeds 
the equipartition field by a factor > 6. Such strong magnetic field 
has been observed only in simulations with 128^ grid points res- 
olution. The rise velocity of the buoyant field in the bulk of the 
convection zone, estimated from z-time "butterfly" diagrams, 
has been found that ui^ ^ O.G^Xrms- From the test-field method 
results, the maximal downwards pumping velocity of magnetic 
field is found to be 7 OAurms- This indicates that buoyancy 
may indeed exist from the base of the convection zone, however, 
the magnetic field expands very quickly during the rise and mag- 
netic structures observed close to the upper boundary are not so 
well organized. 

Besides the emergence of the magnetic field that affects the 
flow pattern locally, other changes are observed due to the back- 
reaction of the magnetic force on the plasma. In the saturated 
phase, the rms-velocity varies between its kinematic value and a 
lower amplitude when the magnetic field is strong. In addition, 
during the peaks and valleys of the field amplitude, the imposed 
shear profile presents systematic variations with respect to its av- 
eraged profile. For most of the models this change is around 3 per 
cent of the shear velocity and is reminiscent of the fluctuations 
observed in the solar rotation profile or "torsional oscillations". 
The variation signal is mainly observed in the places where B is 
strong, with much weaker changes within the convection zone. 
The results encourage the study of the torsional oscillations in 
detail through direct numerical simulations. 
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